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Abstract 

We extensively study the growing behavior of the energy and the pressure compo- 
nents depending on the space-time rapidity in the framework of the Glasma, which 
describes the early-time dynamics in the ultra-relativistic heavy-ion collisions. We 
simulate the Glasma solving the classical equations of motion in the SU(2) Yang- 
Mills theory and systematically investigate the dependence of the Glasma instability 
on the model parameters. We have checked that the transverse and longitudinal grid 
sizes in our simulation are large enough to handle cutoff effects under control. By 
comparing the numerical results from several initial conditions with different mag- 
nitudes of instability seed and also those with different wave-numbers for rapidity 
fluctuations, we clearly see that unstable modes dominantly grow up in the lin- 
ear regime and we also confirm non-linear effects in the time evolution. To extract 
more detailed information on the evolving Glasma, we decompose the energy into 
the components in terms of rapidity wave-numbers. We observe an energy fiow from 
low wave-number modes into higher wave-number modes due to non-linearity in the 
equations of motion. We find that the energy spectrum approaches an asymptotic 
scaling that is consistent with Kolmogorov's power-law form even in the expanding 
system of the Glasma. 

Key words: Color Glass Condensate, Glasma, Relativistic heavy-ion collision. 
Instability, Kolmogorov spectrum 



1 Introduction 



The formation of a quark-gluon plasma (QGP) is achieved in ultra-relativistic 
heavy-ion collision experiments at Relativistic Heavy-Ion Collider (RHIC) in 
BNL and at the Large Hadron Collider (LHC) in CERN. More and more infor- 
mation on properties of strongly interacting quark-gluon matter in quantum 
chromo dynamics (QCD) is being accumulated. However, despite a number of 
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interesting discoveries and new notions established in experiments, the ther- 
mahzation mechanism that leads to the QGP still lacks a deep understanding 
based on QCD first principles. 

There have been many attempts to estimate the thermalization time scale, that 
is phenomenologically known to be Toq < 1 fm/c from analysis in hydrody- 
namic models p!|2] . It was argued [3] that perturbative QCD processes lead to 
a thermalization time of order a~^'^/^ where the saturation momentum Qg 
is a characteristic hard scale in the high-energy reactions; Qs{x ~ y/s) = 1- 
2 GeV in the typical RHIC kinematics. If the prefactor is of order unity, this 
estimate would yield a thermalization time around several fm/c, which is too 
slow to account for the RHIC observations. It has also been suggested that 
QCD perturbation theory, no matter at which order one truncates it, may 
be unable to describe thermalization at all Plasma instabilities, particu- 
larly the Weibel instability [5fS|7|5f^lUfllfl2] . could be promising candidates 
for the early thermalization mechanism [13] and have been well investigated 
in the hard- loop formalism p!^15|16|17|18j . There are also various other ap- 
proaches |T9^,'20"21"22] including a framework based on the AdS/CFT corre- 
spondence [23.24.25.26j. in which thermalization appears as the formation of 
a black hole horizon in the dual theory. 

Parallel to the analysis on plasma instabilities, an important observation was 
made in a related theoretical framework. The occupation number of gluons at 
high energy is so huge that color field amplitudes ~ aj^^'^ can be realized, 
and then non-linearities must be taken into account to all orders. Instead of 
a plain perturbative expansion, then, an approximation in terms of classical 
field equations becomes more suitable because it effectively makes a resumma- 
tion of large- amplitude fields A^. This non-linear regime, whose natural time 
scale is of order ~ 0.2 fm/c at RHIC, is expected to be prevalent during 
the early stages in heavy-ion collisions. The QCD framework known as the 
Color Glass Condensate (CGC) [27f28] rearranges the perturbative expansion 
in terms of classical fields (that correspond to the resummation of infinite sets 
of graphs), and perturbations on top of these classical fields are of higher or- 
der in as- The CGC was first developed in order to deal with the saturation 
effect in the gluon distribution function [2^301131] . Because we are interested 
in the early time (r < Q^^) dynamics before thermalization, the CGC de- 
scription seems to be the most suited for our purpose. The initial condition 
just above the forward light-cone (i.e. at r = O"*") for the classical field that 
appears at lowest order in the CGC description can be nicely formulated in 
the Bjorken (expanding) coordinates [32E3]- Because in the regime of large 
occupation numbers, classical field theory and kinetic theory essentially de- 
scribe the same dynamics [Mf35] . it is conceivable that the CGC can describe 
the isotropization and thermalization processes qualitatively. However, for this 
program to have a chance to work, it is necessary to go beyond the lowest or- 
der description, since at zeroth order it is well known that the pressure tensor 
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never isotropizes. 



It is unfortunately impossible to obtain an analytical solution of the classical 
Yang-Mills equations of motion for the heavy-ion collision (two-source) prob- 
lem, and therefore one has to resort to numerical simulations. In Refs. [36|37f38] 
heavy-ion collisions were simulated in the McLerran-Venugopalan (MV) model 
(i.e. the CGC at lowest order, with a local Gaussian distribution of color 
sources), formulated gauge-invariantly in terms of link variables. These nu- 
merical calculations were later confirmed by an independent simulation in 
Ref . |3n] • The peculiar feature in the CGC description at leading order is that 
the initial condition for the classical fields is independent of the space-time 
rapidity variable r], and the equations of motion maintain //-independence 
during the time evolution. Therefore, the space-time evolution is entirely in- 
dependent of 7], that is, invariant under boosts in the direction of the collision 
axis. Intuitively, one can expect a longitudinally extended structure of the ini- 
tial fields, which is indeed a characteristic property of the so-called "Glasma" 
state [lO], which could be a source for the so-called ridge structure seen in the 
two-particle correlations at RHIC and LHC [¥T]H2] . 

Within the CGC framework, an instability has been found when the boost 
invariant Glasma fields are disturbed [^3)1^ . This instability has some simi- 
larities to the Weibel instability in anisotropic plasmas, but the origin of the 
instability is not yet fully clarified, though several interpretations have been 
proposed P5]|46ll47f48f49j . Also, it is still an issue how to properly formulate 
disturbances to the boost invariant CGC background. However, it is clear that 
such perturbations arise when one considers higher order corrections to the 
leading CGC picture [50] , and therefore they should not be ignored since they 
can alter phenomena such as particle production [5T]|52|53] . etc. Moreover, 
without including these perturbations, one cannot quantify the thermaliza- 
tion processes. 

One of the reasons why the origin of the Glasma instability lacks a clear 
understanding is that the available numerical data from earlier works P3]H^ 
is rather limited. This motivates us to carry out a more systematic survey 
of the Glasma instability, exploiting a wider range of simulation parameters. 
In this work we will specifically investigate the dependence of the Glasma 
instability on the following parameters: 

[i] Transverse grid size N: It is known that the initial energy density at 
r = 0"'" is sensitive to the ultraviolet and infrared regulators in the transverse 
direction, and this sensitivity becomes milder at finite r [1|45II54] . It is therefore 
important to check whether continuum-limit results can be obtained reliably 
by looking systematically at the dependence on the transverse lattice spacing 
a = L/N, where L is the length of the box and the number of lattice points 
in the transverse direction. 



3 



[ii] Longitudinal grid size Nr,: Physical results should also not be contami- 
nated by ultraviolet singularities in the longitudinal direction. To avoid this, 
we will use initial fluctuations with sufficiently small longitudinal wave-number 
V iy is the Fourier conjugate to the rapidity 77). We will conffim that there 
is no dependence on the longitudinal grid size A^,, or the longitudinal lattice 
spacing = L^/A^^, at least until the longitudinal spectrum extends into the 
ultraviolet regions at very late time. 

[iii] Magnitude A of the seeds: In general the strength of the instability 
increases with the magnitude A of the initial ?7-dependent perturbations. This 
dependence is found to be very simple until the unstable modes become large 
enough to be self-interacting: at moderate times, the finite-i^ modes obey the 
linearized dynamics and their amplitudes grow linearly as A increases. 

[iv] Initial wave-number uq: The time evolution depends on the initial con- 
dition in the longitudinal direction, and in particular on the wave-number 
of the seeds. We will show that the instability produced by the smallest wave- 
numbeiQ (z/q = ±1) is the fastest and strongest in a wide range of times. We 
will also conffim the dominance of linearity and the appreciable presence of 
residual non-linear effects by verifying that the result from an initial condi- 
tion with multiple z/q's can be approximated as the superposition of the results 
from initial conditions with single wave- numbers fo's. 

In this paper, we report on numerical results about these dependences for the 
growing behavior at the late stage of the Glasma evolution. These numerical 
results provide a pool of observations that we hope will help to understand the 
microscopic mechanisms that drive the Glasma instabilities. In the section [21 
we recall the most important formulas for a lattice study of this problem. 
In the section [3l we briefly expose some results in the purely boost-invariant 
case. These are well-known results already, and merely serve as a check of 
our implementation. In the section HI we describe what happens when the 
Glasma fields are perturbed by a rapidity-dependent seed, and we investigate 
thoroughly how the result depends on the parameters listed above. We explore 
the time evolution of the i^-spectrum in the section [51 Finally, the section [61 is 
devoted to a summary and outlook. 



2 Lattice formulation of the problem 

In this section, we summarize our equations, with emphasis on the lattice 
specificities. The main goal of this section is to avoid confusions due to the 



^ In our conventions, the wave-number is defined to be an integer between —N^/2 
and +Nn/2. The lowest non-zero wave-number is thus v = ±1. 
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existence of different conventions in the literature for some quantities. 



Tlie link variables are defined as U^{x) = e'^^"'^'^^^^ with the gauge potential 
Afj_{x) that connects the neighboring sites from x and x + fi. This expression 
for the link variable is the most important one in order to establish the corre- 
spondence between the continuum and the lattice (link-variable) formulations. 
Under a gauge transformation by V{x) the link variables transform as 



U,{x) ^ V{x)U, 



and the chromo-electric fields as 



E\x) -> V{x)E\x)V^{x) 



(2) 



The plaquette variable is then defined as 



Uf,y{x) 



U,{x)U,{x + fi)Ul{x + v)Ul{x) 



exp 



X 



(3) 



where the last form is an approximation which gives rise to the field strength 
d/j^Ai, — dpA^ — ig[A^^ A^] of the continuum formulation. This is a very 



useful expression to study the continuum limit. The lattice spacing should be 
the longitudinal one, i.e. a^, which is distinct from the transverse one a, if 
the longitudinal direction is involved in the displacements /}, or i> in the above 
formula. The canonical momenta in terms of link variables are expressed as 



drUi{x) 



-1^ 



E\x)Ui{x) 



drUr,{x) 



-igarjT E'^{x)Ur,{x) 



(4) 



Here we note that in the above we already use variables made dimensionless 
by the transverse lattice spacing a, and thus a will never appear explicitly in 
our formula. Because ?7 is a dimensionless number, does not bring any mass 
dimension unlike a. 

We can also discretize the time in the above equations in order to solve the 
time evolution numerically for Ui{x) and Urf{x); 



U:{r") 
UJt") 



exp 
exp 



-2At ■igE\T')/T'\Ui{T 
-2At ■iga^r'E'^ir'] 



(5) 
(6) 



where t' = t + At and r" = r + 2 A. Note that we have not written explicitly 
the position arguments in the quantities that appear in these equations. Here, 
the exponentiation of E^ is important; otherwise the updated U^ij") would 
not be a unitary matrix. The equations of motion are discretized in the same 
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manner as 



EUt') = E'(t- At) + 2At 



2ga'^T 



Ur,i{x) + U-r,i{x) - (h.c.) 



+ 2Ar^E + f/-,.(x) - (h.c.) 



(7) 



for the transverse components and 



E^ir') = E\t - At) + 2Ar 



2ga„T 



E 



V j=x,y 



Ujr^ix) + U.jr^ix) - (h.c.) 



for the longitudinal component. We recall again that all variables including 
T and At above are dimensionless, and expressed in unit of the transverse 
spacing a. One can make sure that these equations are equivalent to the or- 
dinary equations of motion in the continuum limit, by using the approximate 
form (E]). 

The initial conditions are derived from the requirement of avoiding the singu- 
larity at the collision point. It is written in terms of ul"^\ the classical solution 
of the MV model for a single color source (the label m = 1, 2 indicates which 
of the two nuclei produces the corresponding field). This single-source classical 
solution is a transverse pure-gauge, 



(9) 



where the gauge rotation matrix can be written as 



X 



exp 



1^ 



(10) 



Here, A*^™) is the solution of the Poisson equation, d\A'^'^\xi_) = —p^'^\x±). 
In the MV model, the random color source p^"^^ is Gaussian distributecf^. with 
a two-point correlation 



Higher-point correlation functions are obtained as products of two-point con- 
tractions, via Wick's theorem. In the MV model, p is the only dimensionful 
parameter and it is related to the saturation momentum Qg. Then, the initial 



^ A Gaussian distribution, albeit one with non-local correlations, is also an approx- 
imate solution of the JIMWLK evolution equation [55] . 
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conditions are |36] 



(2)t 



—1 



)t 



+ [UjXx-^x^) - l)[Ur\x-^x. 



t=x,y - 



U^^^Hx-Ax, 



(h.c. 



(12) 



, (13) 



written here in the shghtly simpler case of the SU(2) color group. These dis- 
cretized equations of motion and initial conditions completely define the nu- 
merical MV model. Our conventions are mostly consistent with those of pre- 
vious works, except for the fact that the field is not treated as an adjoint 
scalar-field variable but instead we also describe it in terms of a link variable 
Urj. We note that our formulation is reduced to the conventional one with the 
adjoint scalar-field variable if we take to be sufficiently small. In addition, 
since our interest lies primarily in the qualitative behavior of the Glasma in- 
stability, we disregard the color neutralization that has been imposed onto the 
MV model color sources in other works 

It is now straightforward to implement the numerical calculation for gauge- 
invariant observables such as the components of the energy-momentum tensor. 
They are in the continuum convention written down as 



tr 



L 



_ / rpXX _|_ rpyy 
2^ 



tr 



tr 



Ei 



E^ +B' 

T T 



E' 

L 

Bl 



B' 



E'-B' 



(14) 
(15) 
(16) 



where we define 



^2 ^ ^Va^va ^ 



E' 



LfE'^'^E' 



(17) 



These formulas are implicitly summed over the color index a, for E"^ and 
E^ to be gauge invariant. Here, the chromo-magnetic field squared should be 
expressed in a way consistent with Eq. ([3]), that is. 



xy 



B' 



(ga^r) 



t=x,y 



These definitions are, however, not very convenient when we decompose the 
energy into its Fourier components, as we will see later. We confirm that the 



energy-momentum tensor is traceless; T^n = 



It should also be mentioned that our definition of is different from that in 
Refs. |13pi] by a factor 2 so that Pj, = for an isotropic system. 
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3 Boost Invariant Expansion 



We first check the consistency with the known results in the boost-invariant 
case, i.e. in (2 + 1) dimensions. This is useful also for later discussions on the 
spectral decomposition when we consider ?7-dependent fluctuations. In this 
section, we focus on ?7-independent classical solutions, with initial conditions 
specified in Eqs. (fT2|) and ( IT3|) that are themselves independent of 77. Even- 
tually, however, our goal is to superimpose r^-dependent fluctuations to these 
boost invariant initial conditions in order to incorporate quantum effects [50] . 

We note here that our numerical calculations are limited to the SU(2) color 
group, so that we can get more samples and thus better statistics for the 
same computational cost. As discussed in Ref. [IB], the SU(2) calculations 
capture the essential properties of the Glasma instability, and we expect that 
the results we report here for the SU(2) case would be qualitatively the same 
for SU(3). 



3.1 Initial Configurations 

Let us first take a look at a typical fixed initial configuration, without taking 
the ensemble average over the color sources, in order to develop an intuitive 
picture for the Glasma state. In the MV model, the initial color distribu- 
tion in the transverse plane is akin to white noise, having no correlation be- 
tween different transverse site^. Since the solution of the Poisson equation 

) involves a convolution with a massless (i.e. long 
ranged) transverse propagator, A(a:;_|_) has a smooth structure as shown in the 
left panel of Fig. [H 

Although A{x±) is fairly smooth due to the regularizing effect of this con- 
volution, the corresponding gauge field ^^(a^^), which is a phase of Ui{x±) 
defined in Eq. iQ, is far from smooth. This is actually a consequence of the 
pure-gauge form; contributions from globally smooth parts in A{x^) cancel 
each other and the resulting ai{x±) fluctuates strongly. In other words, ai{xj_) 
has a rough structure because a derivative of A{x±) undoes the regularization 
provided by the convolution. The ai{x±) shown in the figured] is the gauge 
field from a single source only. In a collision, the initial fields such as E^{x±) 



^ This property would be slightly altered if one imposed to the color distribution to 
be color neutral over patches corresponding to the size of a nucleon, or if one were 
using distributions of sources evolved with the JIMWLK equation - in the latter 
case, it has been shown that the JIMWLK evolutions induces correlations among 
the sources over transverse distances of the order of Qj^- 
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Fig. 1. An initial configuration generated by a random color source p{x±). The 
MV-model parameters are chosen as = 64 (number of transverse sites) and 
g^fj.L = 120. Left: The solution A{x±) of the Poisson equation. Only the first color 
component A^{x±) is plotted. Right: Corresponding gauge field ax{x±) calculated 
from e~*^^(^^)e*^'^^'"^^*) = exp[—igai{x±)]. Only the first color component al.{x±) 
is plotted. 

and B^{x±) involve a product of two ai{x±)^s from two independent sources, 
and thus have an even rougher structure. 

It is worth mentioning that the "color-flux tube" structures, whose transverse 
size is expected to be of order Q~^, do not emerge manifestly from the MV 
model since the source distribution has no correlation length in the transverse 
direction. This means that the Nielsen- Olesen instability |48] is unlikely to 
occur at least within the MV model. Nevertheless, we expect that color-fiux 
tube structures would appear after the JIMWLK evolution from the MV- 
model distribution. In the mean-field approximation, in fact, the transverse 
correlation function of color source has been already studied in Ref. [55], which 
should be taken into account in the future works. Also, to make a clear com- 
parison with the Weibel instability and the hard-loop estimates (especially 
early works [T^IITS] ). it would be useful to simplify the MV-model setup in 
(1+1) dimensions without transverse dependences in such a way not to lose 
the qualitative features of the Glasma instability that we will see in the sec- 
tion m In this sense the results reported in the present paper should serve as 
a foundation for future investigations in these directions. 

3.2 Time Evolution 

Figure H] shows the time evolution of the energy density f lT7|) and flTHl) . We 
set g^fiL = 120 which means that g'^fiR^ = 67.7 with nR\ = according 
to the convention for the RHIC physics in the literature [39]. This choice 
corresponds to Ra ~ 7 fm 1.2 x (197)^/^ fm (i.e. the radius of the Au 
atom)] and g"^^ ^ 2 GeV (where g = 2 a.s usually chosen, meaning that 
as — 0.3). In most of the numerical results presented in this work, we keep 
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Fig. 2. Left: Energy density (|14p for various choices of the transverse grid size A^. An 
ensemble average is taken over 30 random configurations of p^"^\x±) for N = 700 
and 300, and over 50 configurations for N = 100, 32, 16 and 8. The MV model 
parameter is chosen as g'^fJ.L = 120. Right: Chromo-electric and chromo-magnetic 
fields ()17p and (jlSp for the same MV model parameters and N = 700. 



0.05 




Fig. 3. Longitudinal pressure defined in Eq. (jl6p and transverse pressure P^, 
defined in Eq. (fT5]l . multiplied by g^^ir. 

using these parameters. We note that the strength of the Glasma instabihty 
as seen in the next section is set by the instability seed A relative to g^fi 
which characterizes the CGC background. In this sense, for our purpose in 
later discussions, the precise value of g^fiL is of no qualitative importance. 

In Fig. 121 we investigate the dependence on the transverse grid size A^: The 
left plot shows how the energy density increases with increasing A^, varying 
from = 8 to = 700 from the bottom curve to the topmost one. It is 
already known in the literature [39] that the A^-dependence (or the UV-cutoff 
dependence) is significant at g'^fir -C 1, but becomes almost irrelevant for 
g'^fiT > 1 (except for small lattice sizes A^ ^ 100). 

As is well-known from earlier works [5B|57] . the energy density decreases as 
r~^, which is the scaling behavior of a free-streaming expansion. In fact the 
energy density multiplied by the time, re, saturates at a time g'^fiT ~ 1 (which 
is sometimes called the formation time, i.e. the time at which partons are 
freed from the wave- function of the nuclei) where the free-streaming expansion 
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starts. This is natural on physical grounds, because the classical fields keep 
interacting strongly until g^fiT ~ 1, and become quasi-free at ^ 1. 

We see that the longitudinal chromo-electric and chromo-magnetic fields be- 
have in the same way, which is the characteristic feature of the Glasma initial 
state [10] • In our convention E"^ , , E"^ , and are all approaching a com- 
mon value for g^fiT > 1, and the system still remains anisotropic (isotropy 
requires that E^ = 2E'^ and B^ = 25^). This becomes clearer if we consider 
the transverse pressure and the longitudinal one defined in Eqs. (ITS!) 
and fll6p . respectively, as shown in Fig. [3J Then, we can find that tP^ and 
tPj^ both approach an asymptotic value and rP^ is much smaller than tP^. 
If Pj^ is exactly zero, this means that the system reaches the limit of the 
free-streaming expansion leading to e oc r~^. If the system is completely ther- 
malized, on the other hand, the pressure must be isotropic and thus P^ = Pl- 
Because an expansion with positive P^ produces work against the expansion 
of the matter, the energy density decreases faster than in the free-streaming 
case, e.g. e oc r~^/^ if P^ = e/3 and if the expansion is purely longitudinal 
(i.e. Bjorken expansion). Therefore, the (leading order) Glasma state alone 
cannot reach isotropization; a candidate for causing the longitudinal pressure 
to grow is the instability due to the ?7-dependent fiuctuations that occur in 
higher order corrections. 



3.3 Spectral Decomposition 



Let us now consider the time evolution of the spectral composition of the 
energy content. That is, the energy is written in terms of Fourier components 

as 



(19) 



(27r)2 ^ J (27r)2 

in terms of the continuum variables. If written in terms of the lattice variables, 
the transverse momenta take discrete values ki = 2imi/L and the integration 
is replaced by a summation over n^. The energy density in the mode k±_ is 
given by 



e^{ki_) 
e^{k±) 



P'"^(-fcx)^'"'(fc±) + T-^(^E'^{-k^)E'''{k^) 
fi''"(-fcx)5''"(fc±) + r^^(P"(-fc±)P*"(fc±) 



(20) 
(21) 



where E'^'^{k±), £'*"(fcj_), B^"-{k±), and B^"'{k±) are the Fourier transforms 
of the chromo-electric and chromo-magnetic fields. The above expressions are 
given in terms of the continuum variables, but it is non-trivial how to re- 
express them in terms of the link variables: Equation ( I2UI) for the chromo- 
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electric component is not changed, but we cannot use Eq. ( 12T|) as it 10 



It is of course possible to introduce the chromo-magnetic field using the rela- 
tion ([3]) in such a way that 



L 



xy 



vy 



-tr 



1^ 



igar, 



-tr 



r 1 - U, 



vy 



pa 

TJX 



-tr 



1 - a, 



rjx 



■ (22) 



This approximation works for sufficiently large number of the transverse sites, 
i.e. > 100. Then, we in principle have two options here: we abandon the 
information on eg{kj_) and just evaluate ej^{kj_) to see the spectral pattern, 
or, we adopt a sufficiently large value of so that we can use the approx- 
imation 022]) • In this work we take the former option, because simulations 
with > 100 would be too time-consuming in the forthcoming studies that 
include //-dependent fiuctuations (and where one thus needs to solve (3 -|- 1) 
dimensional equations). 

At r = 0, the initial spectrum is calculable analytically in the MV model, 
leading to the expression |58] . 



3(^?V)^ 
27r^2 



:ln 



4m2 



k"^: + 4m? + k_i 



/c^ -I- 4m2 — k 



(23) 



with an IR cutoff m. This formula means that ej^{k±) behaves as /c_|_^ ln(A;_|_/m) 
at large k±, which is the expected perturbative tail. 

At non-zero r, we have evaluated ej^{k±) numerically, with results shown in 
Fig. m To draw the ffgure, we have calculated €j^{k± = ^k^ + ky) first as a 
function of k^ = 2TTnx/L and ky = 2'Kny/L and then have taken an average 
over the results that lie in a bin, A^ < y^n^ -|- < A^ + 1, to obtain ej^{ki_) 
with k_[_ = 2txN/ L. We use logarithmic horizontal and vertical axes in the plot 
of Fig. m to facilitate the identification of the initial-time power-law scaling 
in the spectrum. 

We can notice that the energy spectrum at the initial time g'^^r = 0.01 (where 
we started the simulation) obeys the perturbative scaling ~ kj^ up to the UV- 
cutoff, as indicated by a line for eye guide in Fig. H] - in agreement with the 
perturbative expression in Eq. ( 123|) . The spectral shape quickly approaches 



^ This is an issue only because we want to decompose the magnetic energy in 
Fourier modes. If we stay in configuration space, the expressions for and in 
Eq. (jlSp are perfectly fine to calculate the chromo-magnetic energy - the problem 
is that they are not simply a sum of squares. 
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Fig. 4. Energy spectrum Te^{k±) at ^^/ir = 0.01, 0.073, 0.52, and 1500. The trans- 
verse size is chosen as = 300 and the MV model parameter as g^fiL = 120. The 
spectrum at g'^^T = 0.01 shows the perturbative scahng £^{k±) oc /cj^. An ensemble 
average is taken over 30 configurations. 

the asymptotic one within g^fir ~ 1 and it hardly changes after that time. In 
fact the shape at g'^fir = 0.52 is already close to that at g'^fiT = 1500 (where 
we stopped the simulation). Thanks to this, we do not show the results at any 
intermediate time between g'^fir = 0.52 and g'^fj^r = 1500 in order to keep the 
plot legible. 

In later discussions in the section [5] we will perform a similar analysis in the 
longitudinal direction to find non-trivial asymptotic behavior of the spectral 
shapes. 



4 Breaking the Boost Invariance 



After this summary of the boost-invariant results, including the time evolution 
of the energy spectrum, let us proceed with the analysis of the instability with 
respect to //-dependent fluctuations. We will first confirm the existence of the 
instability in a similar setup as the one used in previous studies. For this 
purpose, we choose the initial fluctuations according to the prescription used 
in Refs. [^SlHl] and we fix the MV model parameter g'^fi in the same way as 
in the boost-invariant situation. We take the extent in the t] direction to be 
= 2 units of rapidity and we use periodic boundary conditions for rj as well 
as for transverse coordinates. We note that this extent covers roughly the mid- 
rapidity region where the boost-invariant plateau is expected, and therefore 
it is legitimate to have a boost-invariant background. That is, our parameter 
choice reads; 

2 120 , 2.0 

g fia = — {g = 2) , = — , (24) 

and we perform calculations at various values of the transverse grid size 
and the longitudinal grid size iV^ later. 
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Introducing an arbitrary function /(r/), we can write the ?7-dependent chromo- 
electric fields in such a way that Gauss's law is satisfied by construction, 



6E\x) = a-' [fir] - a,) - /(r/)] CM , (25) 
SE'^ix) = -fiv) E [uRx - - - - Ci^i-)] ■ (26) 



i=x,y 

Before applying Eq. f l2B]) we can also add ?7-dependent fluctuations in f/j. For 
the moment we will discuss the case with SE^ and SE^ only. Here the choice 
of C{^±) is arbitrary and we choose it as a random variable in the transverse 
plane, i.e. 

{C{x^)e{x'^)) = S''S^'\x^-x'^). (27) 

It is important to note that ^* has the dimension of a momentum, which is 
obvious from the above equation. If all the dimensional quantities are scaled 
with the transverse spacing a as is the case in the numerical implementation, 
^* scales as a~^, which is artificial. The strength of the longitudinal disturbance 
should be given independently. To cancel this artificial a dependence, we use 
a dimensionless ^* satisfying {dx±)^-'{x'j_)) = S'^^x,x' (on the lattice) and 
make /(?]) scale as A^~^, so that the fluctuations are proportional overall to 
{Na)~^ = which is a constant. 

In the literature [l3pi] /(?]) is chosen as a random variable too containing 
all infrared and ultraviolet modes in the spectrum. Here, in order to study 
the instability behavior in a well-controlled situation, let us perturb the sys- 
tem with a single //-mode, and consider superposition of multiple modes only 
lateS Hence, we will adopt the following simple form for the moment!^. 

/(^) = Acos(^ry) , (28) 

where = N^ar^ is chosen to be 2 as we explained before and A should scale 
as 

A = ^, (29) 

to make the fluctuations insensitive to the way we discretize the transverse 
plane. We will later fix Aq and then vary N to study the (unphysical) depen- 
dence on the transverse grid size. 



^ It is legitimate to do so as long as the perturbations are small compared to the 
background field, so that their dynamics remains linear. After that, the different 
fluctuation modes will evolve non-linearly and mix. We will return to this point 
later. 

^ In this lattice parameterization of /(?/), the frequency vq is an integer between 
-iV^/2 and -|-iV^/2. 
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4-1 Dependence on the transverse grid size N 

Here we fix uq = 1 in Eq. fl25]) because, as we will confirm later, this lowest 
non-zero mode leads to the fastest and strongest growth in the longitudinal 
pressure. 

Figure O shows how the instability occurs in the longitudinal pressure in 
the long-time run. The left-panel is the pressure multiplied by time [g'^rP^) in 
the unit of {g'^^Y- We note that the plotted quantity is for z/ = z/q = 1 which 
maximizes the Fourier transform of rP^ (except for the zero mode = 0). 
That is, we define, 

P, d'^x j^' dr/ P, (r/, x^) e'(2-/^.). , (30) 

which is complex valued in general^. We thus plot its modulus, |(^^rP^(i^ = 
vq) / [g"^ in the left panel of Fig. O corresponding to an initial disturbance 
at I/O = 1 and Aq = 32 x 10"^ in Eqs. ([28]) and ([29]) (i.e. we choose = 32 
as a reference point and then A = 10"^ for = 32). The horizontal axis 
represents a/^^/tt, which is a natural variable in a longitudinally expanding 
geometry. Note that the Fourier transform P^iv) at finite v is not a measur- 
able quantity in experiments. Its real and imaginary parts both fluctuate from 
negative to positive values depending on the initial conditions and they would 
be vanishing if we compute an average without taking the modulus. Only the 
z/ = mode has a physical meaning as the longitudinal pressure. For our pur- 
pose of studying the instability, we can nevertheless consider {\P^{y)\) (i.e. we 
compute the modulus first, and we perform the average over initial conditions 
next) in order to confirm some known results reported previously [13]H^ . In 
later discussions, we will also decompose the energy density into v components, 
which is a more well-defined quantity to see the instability. 

From the left panel of Fig. the exponential growth with the square root of 
proper time is obvious. The onset of the instability depends on the transverse 
size but the instability slope is rather insensitive to TV. With increasing A^, 
the instability occurs earlier, and eventually the shift in this onset position 
saturates for A^ ^ 32. This means that the instability does not start earlier by 
increasing A^ further. So, we stop increasing A^ and conclude that simulations 
at A^ = 32 are reasonably close to the continuum limit. In the right panel, 
we see how the fields, initially localized in the modes z/q = (background 
field itself) and z/q = ±1 (perturbation), populate the higher z/ modes due 
to the non-linearities in the Yang-Mills equations. For a small perturbation, 
the amplitude in the higher modes decreases very fast (exponentially with z/), 

^ We note that u in our definition ranges from to and the correspond- 

ing wave-number is kri = ^-nvjL^. 
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Fig. 5. Left: Longitudinal "pressure" at = z^Oi where I'o = 1 is chosen for the initial 
condition ()28p . The results are for N = 8, 16, 32, 64, 128 in the transverse direction 
(from the top to the bottom) and = 128 fixed in the longitudinal direction 
(except for = 64 for = 64 and iV,, = 32 for = 128). The seed magnitude 
is A = 10^^ for N = 32, i.e. Aq = 32 x 10^^ in Eq. ([29]). An ensemble average is 
taken over 100 configurations. Right: Spectrum of the longitudinal "pressure" as a 
function of u at g^/ir = 1000 {\/g^J^ ~ 32). 
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Fig. 6. Left: Longitudinal "pressure" at v = uq = 1. The results are for A 



10" 



and A^ = 32 fixed in the transverse direction and A'^ = 128 and A'^ = 256 in the 
longitudinal direction. An ensemble average is taken over 100 configurations. Right: 
Spectrum of the longitudinal "pressure" as a function of v at g^^iT = 1000. 

because in order to produce a field at a given z/, one needs to multiply v seeds 
at vq = 1. 



4-2 Dependence on the longitudinal grid size Nj^ 



Physical results must also not depend on how we discretize the longitudinal 
coordinate. In order to see the dependence on the longitudinal grid size, we 
make a plot in Fig. |6] in the same way as in Fig. [5l where we compare two 
longitudinal grid sizes that differ by a factor 2. As is clear from Fig. [6] we 
cannot find any sizable dependence on N^i, which means that our numerical 
calculations are safely free from lattice artifacts in the longitudinal direction. 
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Fig. 7. Left: Longitudinal "pressure" at v = = 1. The results are for A = 10"^, 
10~^, 10~^ and a fixed lattice size 32 x 128. Right: Spectrum of the longitudinal 
"pressure" as a function of v at g^^T = 1000, -which has a UV cutoff at fmax = 64 
for the A^n = 128 calculation. 



The influence of A'^ can be seen more clearly in the spectrum of |P^(z/)| shown 
in the right panel of Fig. El A different choice of A^^ with fixed = Nr^ar, 
changes the UV cutoff (the spectrum ranges between — A'^/2 and +Nr^/2): in 
the right panel of Fig. [6l v varies from —64 to +64 for the A'^^ = 128 calculation 
and from —128 to +128 for the A^^ = 256 calculation. Since the seeds are at 
V = = 1 m our initial condition, we naturally anticipate that the UV cutoff 
aX, V = Nrf/2 does not affect physical properties around u ^ uq <^ A^r?/2. 

One may think that the results for a larger seed amplitude A = 10~^ in 
Fig. [6] could be sensitive to the choice of A^,, because all the amplitudes in the 
spectrum are substantially larger and thus a change in A^^ leading to a change 
in the UV cutoff may have an influence on the spectrum. We have also checked 
this and found that it is not the case. As is clear from Fig. [6], the difference 
between the results for 32 x 128 with A = 10~^ and that for 32 x 256 with 
A = 10~^ is almost invisible. 



From all these results we can conclude that a value of A'^ = 128 is sufficiently 
large as long as uq is set to be of order unity. If uq takes a larger value, as we 
will see soon below, the spectral shape has a wider distribution and the UV 
cutoff at z/ = A^,;/2 could influence the results. 



4-3 Dependence on the magnitude A of the seed 



In this subsection, we study the dependence on the seed magnitude A of the 
longitudinal fluctuations in Eq. ( l28i) . In Fig. [71 we show some of numerical re- 
sults for different values of A. It is interesting that the results obviously scale 
as I oc A up to the point where the instability growth eventually saturates. 
This A dependence is almost trivial, which suggests that the non-linear effects 
in terms of finite-z/ modes are small unless the unstable modes substantially 
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develop at large g'^HT. One might have thought that any physical observables 
should be expanded in even powers of A and thus |P^(t'o)| should be propor- 
tional rather than A. This is not the case because we took the modulus 
of complex P^{vq) before performing the average over initial conditions, thus 
preventing the linear term in A from vanishing. 

We can also notice from the right panel of Fig. [7] that, once the instability at 
p = i/q gets saturated, the mode aX u = pq stops growing but the spectrum 
spreads quickly to higher z/-modes. Once this happens, the scaling property 
with A is lost. Actually the amplitude of the v = mode reaches one third 
of the zero-mode amplitude (for A = 10~^ results at g'^^r = 1000) before this 
saturation occurs. Then the non-linearities or self-interaction effects cannot 
be neglected any longer. 



4-4 Dependence on the initial seed wave-number uq 

Now we investigate the z/Q-dependence of the instability behavior. First, we 
have chosen three different z/q's a.s uq = 1, 5, and 10 to uncover the general 
trend. Then we study a more general case, that consists in a superposition of 
these initial wave-numbers. 

The figure [8] shows the numerical results. It is clear from the figure that the on- 
set of the instability is delayed for larger uq. This is naturally understandable 
from the kinetic term for finite-i^ modes. The (positive) kinetic term generally 
tends to stabilize the system, and so the (negative) term responsible for the 
instability is overwhelmed if u is large enough. For dimensional reasons, how- 
ever, this kinetic term for the finite- z/ modes is of the form i^^/r^, and therefore 
it decreases with time. In contrast, the negative term that generates the insta- 
bility has a constant coefficient of order unity - therefore, it dominates over 
the kinetic term at sufficiently large time. In other words, the kinetic term in 
z/^/r^ only delays the start of the instability by a time proportional to z/. We 
can confirm qualitatively this interpretation from Fig. [HI The growth appears 
to start around y/g'^jrr ~ 9.5,21,30 for z/q = 1,5,10, respectively, and the 
starting times are indeed ordered proportionally to uq. 

The spectral shape shown in the right panel of Fig. |8] looks discontinuous for 
z/q = 5 and 10, but this is a trivial artifact of the fact that our initial seed 
contains a single z/q mode. Indeed, the non-linear couplings in the Yang-Mills 
equations can only produce linear combinations with integer coefficients of 
the frequencies present in the initial condition. For instance, with = 5, only 
modes with u a multiple of 5 can be produced by the non-linear evolution. 
The fact that the amplitude if these higher harmonics is suppressed compared 
to the amplitude of the base frequency uq is a good indication that we are still 
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Fig. 8. Left: Longitudinal "pressure" at v = vq with z/q = 1, 5, and 10 in the initial 
condition (I28p . The lattice size 32 x 128 is fixed and the seed amplitude is A = 10~^. 
Right: Spectrum of the longitudinal "pressure" as a function of ly at g'^fiT = 1000. 



Non-linear Result 




-30 -20 -10 10 20 30 



Fig. 9. Left: Longitudinal "pressure" at = 1, 5, and 10 with the initial condition 
given as a superposition (j3T]) . Right: Spectrum of the longitudinal "pressure" as 
a function of u at g'^fiT = 1000. The solid line represents the results from the 
non-linear evolution with the initial condition given as a superposition, while the 
dotted one is a simple linear superposition of the results shown in Fig. [8l 

in a regime where the non-hnearities are rather small. 

We can also see that the non-linear effects are certainly present but not very 
large by repeating the above calculation with an initial condition given as a 
superposition of three modes, namely, 

/(ry) = AEcos(^ry) (31) 

with z/q^^ = 1, ul^^ = 5, Uq^^ = 10, instead of the single-mode fluctuation of 
Eq. fl28l) with uq = 1, z/q = 5, uq = 10 individually. If the non-linearity is 
significant, the results should not be a simple superposition of the individual 
calculation. In view of Fig. [HI however, we can see only small deviations from 
the individual results shown in Fig. |8]at g'^fiT = 1000. This is an important 
observation because the initial fluctuations in reality should have some con- 
tinuous spectral distribution (see Ref. [50] for instance) and we simplified our 
analysis by considering the special case of a single-mode fluctuation. As long 
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as the non-linearity is a minor effect, we can just add up the single-mode 
outputs in order to get approximate results for general initial spectra. 

Although they are several orders of magnitude smaller at g'^^r = 1000, we see 
that higher harmonics mixing the three z/q's are enhanced by non-linear effects 
in the right panel of Fig. [91 For instance, the amplitude of the mode u = 21 
is much larger than the sum of the amplitudes at z/ = 21 for the three seeds 
considered in isolation. Indeed, when the seed z/q = 1 is alone for example, the 
mode z/ = 21 is produced only after 20 interactions, which explains why it is 
totally negligible. Furthermore, the seeds z/q = 5 or 10 can never lead to z^ = 21 
by themselves. When the three seeds are superposed in the initial condition, 
21 can be reached more efficiently as 10 + 10 + 1, i.e. with 2 interactions only. 

Apart from small non-linear effects, it is certainly true that the smallness of 
the field fluctuations we have considered gives a justification for our simplified, 
single mode, analysis. However, at the same time, thermalization is hardly ex- 
pected to occur in these circumstances, precisely because in this regime the 
interaction effects are small. In realistic conditions, the parameter A that con- 
trols the fluctuations is related to the strong coupling as, and therefore it could 
be larger than Aq = 10~^ that we have considered here. As a consequence, we 
expect that the non-linearities would start much sooner and be much more 
important in applications of this framework to actual heavy-ion collisions. 
The study with such large A is beyond our current scope because a large 
A requires full information on quantum fluctuations and needs appropriate 
renormalization procedures, which is technically involved |59] . 



5 Time Evolution of the Longitudinal Spectrum 

In spite of the instability the zero-mode is still dominant among all the modes 
in i\(z/). With an initial condition that contains a single mode as in Eq. (128|) . 
the V = vq = 1 mode is the second largest, which is obvious from the spectrum 
shown in the right panel of Fig. [5l Moreover, the spectrum is spreading into 
higher wave-numbers as the time elapses. In this section, we shall discuss the 
scaling property of the energy spectrum with respect to the frequency z/, and 
show that it possibly displays also a Kolmogorov's turbulent spectrum. A tur- 
bulent spectrum has already been reported in the context of the QCD plasma 
instability piTS] , but these results with the power ^ —2 instead of —5/3 are 
in conflict with Kolmogorov's wave turbulence. On the other hand, the study 
performed in Ref. [SO] has shown that a Kolmogorov spectrum appears due 
to instabilities in Yang-Mills theory, for a fixed-volume system. Therefore, it 
is still worth addressing the turbulent spectrum in the present Glasma sim- 
ulation, that has longitudinal expansion. This is because the Glasma system 
is expanding along the beam axis, which tends to tame the turbulences and 
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Fig. 10. Time evolution of the chromo-electric energy density as a function of the 
scaled longitudinal wave-number v j {g^ ^t) . The lattice size is 32 x 128 and the seed is 
A = 10~^. Left: Results at early times up to g^piT ~ 1500. Note that the horizontal 
axis is not logarithmic and also that the vertical axis is for not re^ (u) but r^e^ (i/) 
(see the text for details). Right: Results at later times up to g^fiT = 3000. 

to delay the isotropization or thermalization in general. It is thus non-trivial 
whether the universal scaling still holds or not in the expanding system. Be- 
sides, our description of the Glasma evolution is purely classical and the precise 
relationship between the Glasma instability and the QCD plasma instability 
is not fully understood yet. 

Two plots in Fig. show the longitudinal energy spectrum (for the chromo- 
electric part only), that is, the energy density decomposed in terms of the 
longitudinal wave-number u as 

f , ^ , , 

^E=hr:eM (32) 



with 



e^[u) = (tr 



27r 



^'"^(-z/)E''"(z/) + r-2(E*"(-z/)E^"(z/)) ). (33) 



(The transverse momentum k± of the fields is integrated over.) The question is 
then what kind of scaling spectrum we can expect for Ej^^u) at later times. For 
this purpose, we consider three characteristic quantities; the wave-number, the 
Fourier decomposed energy, and the rate of energy flow in z/-space (denoted 
by ip below). The wave- number i/ itself is, however, dimensionless and it seems 
difficult to cope with the dimensional matching argument by Kolmogorov. We 
propose that one should use a scaled variable v/{ct) instead of in the case 
of an expanding geometry. This is because err/ is approximately equal to the 
longitudinal coordinate 2; if 77 is small enough, and vjicr) is approximately 
equal to the momentum "p^. Then, the dimension of the scaled wave- number is 
[i//cr] = We note that the integration variable in the decomposition (1321) 
should also be vjcT for this scaling analysis, which means that the relevant 
energy density per mode should be CTSj^iy). This in addition should be mul- 
tiplied by the expanding length cr. Therefore, the dimensions of the three 
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Fig. 11. Time evolution of the chromo-electric energy density as a function of the 
scaled longitudinal wave-number u/{g'^iJLT) up to g^^iT = 4500. Left: Results in- 
cluding later time spectra. The horizontal axis is logarithmic to make it evident 
how the power-law spectrum spreads. Right: Longitudinal "pressure" with three 
regimes identified schematically. 

relevant quantities are 

[u/ct] = l-\ [V^icrfeM] = ^'t'^ M = I't'^ (34) 

which are identical to the standard analysis of Kolmogorov's spectrum. Pro- 
vided that Vj_(cr)^eg(z/) is expressed in a form of {p / ct)'^ {ipY , we can fix a 
and /3 uniquely from the dimensions and thus expect the Kolmogorov scal- 
ing, that is, r^£:^(i/) oc (z//r)~^/^. In fact we can confirm this scaling property 
asymptotically from the right panel of Fig. [101 

One might be wondering what the fate of the power-law spectrum should be at 
even later times. We understand from Fig. [TT]that the Glasma time evolution 
in the present MV-model simulation has three stages; the linear regime, the 
non-linear regime, and the UV-cutoff regime. In the linear regime the unstable 
modes grow up and the non-linear effects are still minor. When the amplitudes 
of z/ 7^ modes become so large that the non-linearities can affect the time 
evolution, the instability stops, as shown in the right panel of Fig. [TTl and 
the non-Abelianization occurs [7fH] . Because non-linearities are no longer 
negligible in this regime, the energy cascade due to non-linearities is naturally 
expected, which would lead to the Kolmogorov-type spectrum. Eventually, as 
the turbulence decays into higher UV modes, the UV-cutoff effect flattens the 
whole spectrum as plotted in Fig. [Tl]up to g^jJiT = 4500. 

To the best of our knowledge the present work is the first to exhibit a Kol- 
mogorov spectrum from the Glasma time evolution. Our results imply that 
the framework of the Glasma, in terms of classical variables together with rj- 
dependent fluctuations, may encompass a rich physics contents including the 
QCD plasma instability. 

Here, as a final remark in this section, we draw attention to related works 
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in analogous systems. In the context of an ultra-cold gas, the Kolmogorov 
spectrum has been obtained as a solution of the Gross-Pitaevskii equation 
(GPE), which is a classical equation of motion like the Glasma equations 
The calculational procedure by means of the GPE is similar to the Glasma 
simulation, but in that case the microscopic mechanism of the energy flow is 
understood in terms of the vortex dynamics. In future studies on the Glasma 
dynamics, it will be important to clarify the structure of the energy flow, and 
it would be very interesting to see whether the so-called non-thermal flxed 
point [52] appears in the strongly correlated Yang-Mills theory and Glasma. 
Indeed, a possible deviation in the power could be expected in Yang-Mills 
theory by non-perturbative resummations [63] . 



6 Summary and Outlook 

In this paper we investigated the time evolution of the classical equations of 
motion in the SU(2) pure Yang-Mills theory, especially in the framework of 
the Glasma or the McLerran-Venugopalan model in (3-1-1) dimensions. Our 
emphasis has been put on the systematic study of the instability behavior 
that occurs in fluctuations with respect to the space-time rapidity r] when the 
initial conditions break boost invariance. 

Computing the longitudinal pressure component at non-zero wave-number z/, 
that is a conjugate of rj, we have numerically conflrmed the Glasma instability 
that had been found in preceding works [13]H^ . We carefully checked the sys- 
tem size dependence of the instability to conclude that the Glasma instability 
is a robust consequence from ?7-dependent disturbances and is not sensitive 
to how the transverse and longitudinal coordinates are discretized. Hence, the 
Glasma instability is certainly a physically well-deflned phenomenon. At the 
same time, however, this means that taking the continuum limit does not help 
with an earlier onset of the instability. The typical time scales we found for 
the Glasma instability is of order of hundreds in terms of g^fiT, which is ar- 
guably too slow for being useful in heavy ion collisions. One should however 
keep in mind that the seeds we have used in this numerical study are very 
small compared to the Glasma flelds, while in a realistic situation they are 
suppressed at most by a power of the strong coupling constant. There may 
also be a chance that the non-linear effects could more or less accelerate the 
instability. We compared two results for the spectral distribution, one from 
the superposition of individual calculations with the initial condition that is 
given by a different wave- number uq, and the other from a single calculation 
with the initial condition that contains all i^o's. We could observe the non- 
linear effects in our numerical results, which would not affect the instability 
behavior substantially, however, as long as the instability seed A is reasonably 
small. 
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We have decomposed the energy density into Fourier modes to see how the low 
wave-number modes cascade toward higher wave-number modes due to non- 
hnearity in the equations of motion. Our motivation to do so is to examine 
whether the so-called Kolmogorov's scaling can appear in the turbulent system 
described by the Yang-Mills theory put in an expanding geometry. In contrast 
to the boost-invariant Kolmogorov-type scaling behavior in the energy 

spectrum a function of the longitudinal wave-number u may appear when the 
boost invariance is broken by the initial fluctuations. In this case, the violation 
of boost invariance increases with time due to instability, but the u = mode 
remains very large compared to the other u ^ modes, which enables us to 
regard the z/ = mode as the source for the energy injection. Then, we could 
clearly observe that the energy spectra approach an asymptotic form that 
shows the Kolmogorov's scaling with the power —5/3 ~ —1.67. As far as we 
know, our work is the first example that exhibits the Kolmogorov-type scaling 
law in the expanding systems. Although it is difficult to find any indication 
that the system indeed becomes isotropic and thermalized in our simulation, 
the confirmation of the Kolmogorov-type spectra is a promising signal for the 
general tendency heading for thermalization. 

There remain very interesting questions to be investigated in the future. In 
order to unveil the microscopic mechanism for the instability, it would be very 
useful to make approximations on the classical equations of motion in such 
a way that the analytical treatment could be feasible. For instance, as we 
have confirmed, the instability mostly lies in the linear regime as long as its 
magnitude specified by A is small enough. Therefore the linearization of the 
equations of motion in terms of fiuctuating fields on top of the boost-invariant 
CGC background should correctly describe the instability in early times. In 
this work we did not exploit such a comparison between the numerical output 
and the analytical result from the linearized equations. This is because it is not 
quite efficient to linearize the equations of motion given in terms of the lattice 
variables. The linearization is of course possible, but it would eventually be 
reduced to similar calculations in the continuum formulation. Then, it would 
make much more sense to make a comparison between the numerical and 
analytical results that are both formulated in the continuum variables. We 
have already executed test programs to make sure if the numerical results 
from the Glasma simulations are consistent with each other, given they are 
written in terms of lattice and continuum variables. We will report on this 
analysis elsewhere. 

It is also an intriguing question to look for a link between Kolmogorov's tur- 
bulent spectrum and the chaotic behavior generally seen in the solution of the 
classical Yang-Mills equations. This problem may have some relevance for a 
deeper understanding of the microscopic description of the entropy generating 
process and how the thermalization is achieved at all. The time evolution of 
the Glasma and its physics contents deserve much more investigations in the 
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future. 
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